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L INTRODUCTION 



A. ELECTRONIC INTELLIGENCE 

Radar has been called the greatest advance in the sensing of remote objects since 
the invention of the telescope in 1608 [Ref. 1], Despite the numerous advantages of using 
radar for locating and classifying remote objects, there is a serious flaw. A large amount 
of energy must be transmitted in order to produce a radar return suitable for processing. 
Electronic Intelligence (ELINT) is the process of observing this energy passively and 
using it for locating or classifying a remote, radio frequency (RF) emitting object. 

Because a radar signal must travel twice the distance from the object to the radar 
transmitter, ELINT processing should carry a range advantage of at least two. Strict 
numerical treatment of the range advantage of ELINT over radar is given in Appendix A. 
ELINT also has the advantage of not requiring any energy emission to accomplish its 
task. This allows the ELINT receiver to remain covert. 

It is relatively easy to intercept a radar signal (they typically are of a high power 
density). If the signal can be detected, then comparison with another receiver will give 
the angle of arrival (AOA) of the signal relative to the sensor pair. This is useful for 
establishing the position of an emitter if it is stationary. If the emitter is moving, further 
processing of the angles of arrival can be performed to estimate the emitter track. 

Despite current trends in the field of Low Probability-of-Intercept (LPI) radars, there 
can be no doubt that ELINT will still be used. The concept of an LPI radar has been 
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described as trying to illuminate an incoming aircraft with a search light without having 
the aircrew able to see the spotlight [Ref. 2], Though LPI radars may require special 
techniques and procedures to be detected, it is unlikely that they will remain immune to 
detection for long. 

B. GOALS FOR THIS RESEARCH 

A simple, time domain algorithm which may be used to process radar data from a 
dual sensor array is desired. The method must extract useful times of signal arrival in 
order to provide a tracking algorithm with angles of arrival of the signal. The method 
must also produce good estimates of primary radar signal parameters. Of specific interest 
for this work are the radar pulse repetition frequency (PRF), the scan rate, and the signal- 
to-noise ration at the receiver (SNR). The precise meaning of the terms PRF and scan 
rate, as they apply to this work, are discussed in the next chapter. Some a priori 
knowledge of the radar pulse width is assumed. 

The method advocated here requires a minimum amount of data to be passed 
between the receiver platforms. This lowers the required bandwidth of the 
communications channel between the receivers and also eliminates the need for a central 
processing station on the ground. Times of arrival are generated by both platforms and 
synchronized against a highly accurate master clock. Times of arrival, synchronization 
information, and receiver position information are all that must be transferred. 
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II. BACKGROUND INFORMATION 



A. BASIC RADAR PARAMETERS 

There are numerous parameters which can be used to describe the characteristics of 
a given radar system. In this section, we will review those of interest for this work. For 
further information regarding radar parameters of interest to ELINT processing consult 
Electronic Intelligence: The Analysis of Radar Signals by Richard G. Wiley [Ref. 1], 

1. Scan Rate/Scan Pattern 

The beam from a radar antenna typically covers only a small portion of the 
angular area of interest. Therefore, the beam is moved over time. Many different 
patterns may be used to sweep out a volume. The scanning pattern is designed to cover 
the volume of interest of the radar. The size of the volume of interest may change over 
time. It is not unusual for a radar sensor to search through a hemisphere in order to 
acquire targets. Once the target has been acquired, smaller volumes will be analyzed to 
perform the function of tracking. The change in scan rate or pattern associated with the 
switch from search mode to track mode or to a separate tracking radar is of interest. 

Another use of scan pattern information is in determining possible capabilities 
of the radar emitter. Radars with mechanically produced scan patterns tend to be 
periodic. Electronically switched scan patterns may or may not be periodic and typically 
can change scan patterns rapidly. Radars with more complicated scan patterns tend to 
have more capability for target localization and tracking than those with simple ones. 
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There are fan beam and pencil beam radars. Pencil beam radars have very 
narrow antenna beams and use complex scan patterns to cover their volumes of interest. 
Typical search patterns for pencil beam radars are raster, helical, and spiral scans. 

The second classification, fan beams, have radar antenna beams which are 
characteristically much wider in one direction than another. This means that the radar has 
very broad coverage in one direction, but poor resolution. Therefore, the radar is moved 
perpendicularly to the plane along the broadest coverage line. A common configuration 
for a long-range, early warning type radar is to have a fan beam radar which scans in 
azimuth. Figure 1 after [Ref. 1] shows the geometry involved as well as the signature 
seen at the ELINT receiver. This arrangement gives good long range detection capacity 
with high accuracy in azimuth, but poor elevation information. A second height- finding 
radar can be used to gather elevation information [Ref 1], This work will concentrate 
on detecting signals from these types of radar installations. These often provide long 
range detection or early warning capability and are of interest. For the remainder of this 
work, we will be interested solely in these types of acquisition radars. The term PRF will 
be used to mean the frequency at which pulses are emitted from the radar. The term scan 
rate will apply to the number of revolutions in azimuth that the radar makes per unit time. 

The signal bursts in Figure 1 are free of background noise, but show the 
essential pattern to be sought in sampled data from an ELINT receiver. Here we present 
a possible method for estimating the energy centroid of such signal bursts in order to 
provide accurate times of pulse burst arrivals. Correlating information from two sensors 
about these times of arrival will give the direction of the signal (DOA). An algorithm is 
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Pulse Amplitude 




(a) CIRCULAR SCAN REGION 




Figure 1. Typical Configuration of a Long Range Search Radar 
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proposed to obtain accurate centroids simply using a decision tree for detecting 
individual pulses inside the burst and then smoothing them with a Kalman Filter. 

The other commonly used pattern for a fan beam radar is sector scanning. 
This is used when only a portion of the azimuth circle is required for coverage. Height 
finding radars generally use a sector scan in elevation when coupled with early warning 
radars. Sector scan radars produce pulse bursts identical to those of circular scan radars. 
However, the bursts will not necessarily come at even intervals. This work will not 
attempt to identify sector scan patterns. 

2. PRI/PRF Analysis 

The pulse repetition interval (PRI) is an important radar parameter to 
evaluate. For a radar with a fixed PRI, the maximum unambiguous range and velocity 
are given by the following equations: 

R u = C X 2 RI ’ (1) 

and 

V = - ; (2) 

u (2 x PRI xRF) 

where 

R^ - Maximum unambiguous range (m), 

V u = Maximum unambiguous radial velocity (m/s), 

PRI = Pulse Repetition Interval (s), 

RF = Radar carrier frequency (Hz), and 
c = Speed of light. 
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These two equations provide more than just information about the emitter platform 
capabilities. Information such as the maximum range of the radar could serve as the basis 
for initial estimates of emitter position For high PRF radars, the unambiguous range is 
very short. Ranging, then, is performed by other techniques such as FM ranging (chirp). 
In this case, Equations (1) and (2) are not valid. 

There are radars with constant PRI and jittered PRI. In this simulation only 
signals of constant PRI were used. The method should be able to estimate an average of 
the radar PRI for all types of radars, including those with PRI jitter. However, the 
resolution of the techniques in this work is not sufficient to accurately measure the 
magnitude of PRI jitter. 

3. Pulse Width/Pulse Shape Analysis 

Finally, we will briefly discuss the importance of pulse shape. For any digital 
sampling method to detect a pulse stream from a radar, the sampling frequency must be 
at a minimum of one sample is taken over the pulse width. If this is not true, then it is 
possible for pulses to go undetected. This will be seen to adversely affect time difference 
of arrival (TDOA) methods. 

B. THE TDOA GEOMETRY 

This section introduces the concepts and equations necessary to convert TDOA data 
into useful estimates of the angle of arrival of a signal. 
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Hyperbolic Lines of Position 



For the TDOA problem, the essential parameter is the arrival time difference 
At of a signal between two (or more) receivers. The resulting time difference measured 
at the receivers could be the result of an emitter anywhere on a hyperbola constituting 
possible positions. This is called the hyperbolic line-of-position (LOP) of the emitter. 
The locus of points on which the emitter could lie for a given At is given in Figure 2. 
For this discussion, both receivers will always be assumed to lie on the x-axis. This 
assumption allows the hyperbola to be easily evaluated and can be obtained by a simple 
coordinate transformation. 

A hyperbola is the set of points in a plane whose distances from two fixed 
points in the plane have a constant difference. The two fixed points are the foci of the 
hyperbola [Ref. 3], For the TDOA problem, the foci of the hyperbola are the receivers. 
The hyperbola can be constructed knowing only the absolute value of the At. If we know 
the sign as well, i.e., which receiver was first to intercept the signal, then one side of the 
hyperbola may be ignored. Letting d be the distance between the receivers the following 
relationship is true: 

At max — / (3) 



where 

d = is the receiver separation distance. 

At max = the maximum possible At based on receiver separation, and 
c = the speed of light. 



8 




Figure 2. The Hyperbolic Line of Position Geometry 
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Once the At max is established, the line of position (LOP) calculation follows easily. The 
center to foci distance is simply d/2, so the center to intercept distance is 



_ _ d At 

<3 — — • 



2 A 



The equation of any hyperbola is given by 



x 



+ yl = 1 ; 
2 b 2 



where 



(4) 



(5) 



b = 




( 6 ) 



Equation (5) is the hyperbola of position for the emitter for the TDOA problem. Figure 
2 shows the quantities graphically. The hyperbolic line of position, as well as its 
asymptotes are shown. It is not evident from these equations, but if the sign of TDOA 
is known, then one curve of the hyperbola can be neglected. From Figure 2 it is clear 
that, at long range, we can approximate the emitter LOP just using the asymptotes to the 
hyperbola. The equation of the asymptotes is 

y = ±-x. (7) 

a 



2. The Angle of Arrival Approximation 

The usual ELINT situation involves long detection ranges relative to the 
receiver separation. Due to this, the correction for hyperbolic lines of position is not 
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necessary. We can use the approximate angle of arrival in order to describe the emitter 
location. Figure 3 shows the geometry after some simplifying assumptions. If the emitter 
is very distant, the approaching wave fronts of the radar signal are nearly straight, parallel 
lines. We can measure the At between receivers. Knowing the distance between receivers 
gives At max as before: 



a = the angle of arrival (see Figure 3). 

This approximation will be used for all LOP calculations. No treatment is given for the 
case of emitters which are too close for this approximation. 




( 8 ) 



where 



11 




Figure 3. The Angle of Arrival Geometry 
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III. RADAR DATA SETS AND THEIR SPECTRA 



A. OBTAINING THE DIGITAL DATA 

In order to use a digital computer to process signal data, the data must be sampled 
to put it in a form which the computer can manipulate. This technology is readily 
available. The special consideration for radar signals is that they are of high frequency 
and, depending on which of the signal characteristics are to be examined, the sample rate 
may need to be very high. Measuring details about the individual pulses, for example, 
would require a sampling frequency greater than the Nyquist rate for the carrier frequency 
of the radar. This would severely tax the data handling capabilities of most computer 
systems. The method used here requires a sampling frequency which provides at least 
two samples per pulse. This turns out to be 

SR z — ; (9) 

pw 

where 

SR = the sample rate, and 

pw = the pulse width of the radar emitter. 

The radar signals of interest are from pulsed type radars. Continuous wave radars 

do not lend themselves to analysis by TDOA methods [Ref. 4], The continuous wave 

radars will show a phase shift at two receivers at different positions. However, if the 

phase shift is greater than 2 tc, the number of phase shifts cannot be unambiguously 
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determined. It is possible for the same phenomenon to affect pulsed radars as well. One 
receiver could be far enough away that it is lagging by more than a single scan cycle. 
This is not likely since typical scan rates are on the order of 1 cycle per second or less. 
For one receiver to lag another by an entire scan cycle, assuming a rate of 1 scan per 
second, would require a separation distance of 300,000 kilometers. 

The typical shape of a pulsed radar with no frequency modulation is seen in top 
portion of Figure 4a. The signal must be envelope detected so that digital sampling can 
occur. The frequency of the carrier wave is in the GHz range and cannot be sampled at 
a Nyquist rate, hence the envelope detection method must be used. By this method we 
now have approximately a square wave with a fixed duty cycle. This cans easily be 
converted to digital data as illustrated in Figure 4a. 

The amplitude of the square wave increases as the radar increases in strength during 
the scan. Peak occurs at the boresight of the radar. One half period of a sine function 
was used in constructing simulated data sets. Figure 4b provides another view of a data 
set with the characteristic features labelled. 

B. TYPICAL SIGNALS OF INTEREST 

The previous section describes the basic procedure for obtaining the digital data 
used. Here we will examine the data more closely. All emitters are simulated using the 
MATLAB programming language. Parameters such as PRF and scan rate were chosen 
by the author to be convenient numbers. Figure 5 shows one burst from a simulated data 
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Pulse Amplitude 



Original Signal 



Envelope Detection 



Sampled Data 

(a) Transition to Sampled Data 




(b) Signal Terminology 



Figure 4. Signal to Data Conversion and Terminology 



set used in this research. This data simulated a radar with a PRF of 9.0 kHz and a pulse 
width of 5 ps. The pulse burst is shown without the presence of noise. The top graph 
shows the entire pulse burst. A cosine function is used to produce the burst envelope (the 
change in amplitude of the received signal as the main beam sweeps past the receiver). 
The second graph shows a closeup of some of the pulses inside the burst. 

We are interested in the frequency components of the typical pulse burst so we will 
examine its frequency spectrum. The power spectral density was computed and is 
displayed in the upper graph of Figure 6. From the first graph we can see that the first 
null in the spectrum occurs at 200 kHz which corresponds to l/(pulse width). In the 
lower graph we observe that the first spectrum line occurs at the pulse repetition 
frequency of 9.0 kHz. These two parameters are readily evaluated from the spectrum of 
the data. 

Figures 7 and 8 are included to show the effect of sampling rate on the pulse burst 
and the raw data set. The plots are produced in exactly the same method as before, but 
at a sample rate of 400 kHz. In the previous example, the sampling rate was 1 MHz. 
As we can see in Figure 7, there are only two samples inside each radar pulse. In Figure 
8 the spectrum shows us that the null occurs at the very edge of the spectrum and cannot 
be used to estimate the pulse width. The methods discussed in the following sections of 
this chapter assume a sampling rate at or very close to twice the pulse width. This 
ensures good detection performance with a minimum of data. This also accounts for the 
statement in the introduction that a priori knowledge of the pulse whdth is required. 
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Figure 5. Pulse Burst (Sample Rate — 1 MHz) 
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Magnitude • Magnitude 



Pulse Burst Power Spectral Density 




Closeup of Pulse Burst Power Spectral Density 




Figure 6. Pulse Burst Spectrum (Sample Rate = 1 MHz) 
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Typical Pulse Burst 





Figure 7. Pulse Burst (Sample Rate = 400 kHz) 
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Magnitude Magnitude 



Pulse Burst Power Spectral Density 




Closeup of Pulse Burst Power Spectral Density 




Figure 8. Pulse Burst Spectrum (Sample Rate = 400 kHz) 
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IV. TIME DOMAIN BASED SIGNAL ANALYSIS 



One of our goals was to produce a method for detecting and locating medium to 
high PRF signals. The method developed uses a decision algorithm to process raw, time 
domain radar data. Estimates of the values of PRF and scan rate are smoothed using 
Kalman filtering techniques. 

The Fast Fourier Transform is used during the initial phase of acquiring a emitter 
in order to produce estimates of the received SNR as well as providing a check on the 
time domain measurement of PRF. The remainder of the processing is carried out in the 
time domain. Optimal estimation theory (optimal filtering of signals embedded in noise) 
is based on time domain methods, as are many aspects of modem control theory (i.e., 
optimal control and adaptive control [Ref 5.]). 

All of the source code files for these algorithms are given in Appendix B. No 
discussion of the coding of the algorithms is presented in this text other than comments 
included in the source code files themselves. 

A. PRF/PRI MEASUREMENT 

The algorithm used to process the radar data depends on a floating threshold. 
Ideally, we would discuss this first. However, it will be much easier to assume that the 
threshold value is known a priori and develop the measurement techniques. Once it has 
been shown how the signal parameters are obtained, the actual method for setting the 
threshold will be discussed. The PRF/PRI measurement forms the central process in the 



21 



overall method. Every incoming bit in the data stream must be processed to find out if 
a pulse has occurred. If it has, then a flag is set and the time between pulses is retained. 
This forms the basis of an estimate of the PRF of the radar emitter. Further logic allows 
us to measure the scan rate and to produce an estimate of the burst envelope to be used 
for the TOA problem. A first order Kalman estimator is employed to smooth out 
individual estimates of the signal PRF. 

As the pulse stream enters the processor, each sample is checked to determine if it 
exceeds a threshold value. When two samples in succession exceed threshold, then the 
algorithm declares that a pulse has occurred. The time between detected pulses serves 
as an observation of the PRI. Since the PRI is assumed to be constant for a specific 
emitter, the following equations serve as the system model for the PRI [Ref. 6]: 



where 



**♦1 = x k ; 

**♦ i = + 



( 10 ) 



x k = PRI estimates, 

z k = PRI observations, and 

v k = observation noise with variance r 0 . 

The discrete optimal filter for this system is 



**♦1 = x k + 



1 lr 0 



* + l/r 0 



(z k - x k ). 



( 11 ) 



If the noise v k is Gaussian noise with zero mean and r 0 =l the optimal filter is 
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**♦. = ** + -^77 (h ~ **)• 



( 12 ) 



In the simulation, we have assumed that the sampled data is corrupted by Gaussian white 
noise of zero mean and unity variance. The noise out of the envelope detector would, of 
course, be Rayleigh distributed. Gaussian noise is used as a simplifying assumption. The 
addition of Gaussian noise to the signal data can also be interpreted as modelling noise 
sources other than the background noise (i.e., thermal noise in the receiver past the 
envelope detector). 

It was noted during development that the noise rejection of the method was not 
adequate. Once the algorithm detects the start of a pulse burst, it changes its operation 
to begin estimating the parameters of interest. After the first detected pulse, the next 
several will be susceptible to being 'dropped.' Dropped pulses are those pulses which are 
not detected and are inside the region of the pulse burst that the method integrates to 
determine energy centroid for use in the TOA problem. 

There are several mechanisms which can lead to dropped pulses. Due to the 
presence of noise one of the samples inside a pulse might be below threshold. It is 
possible to obtain only one sample during a pulse if the sample frequency is too low or 
if the radar pulse is short. Radars are susceptible to skipping pulses. Therefore, there 
may not have been a pulse from the emitter. Nothing can be done about pulses that the 
emitter never sends. If these are infrequent events, the Kalman estimation will still be able 
to produce good results. 
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If a pulse is dropped, the value of the residual in the above equation becomes large 
enough that the PRI observations will be in error. Remember that a dropped pulse yields 
an observation with a relative error of 100 percent. Since there are relatively few pulses 
per burst, (typically from 90 to 150 pulses exceeding threshold), a single dropped pulse 
per burst becomes a significant error. 

To eliminate dropped pulses due to noise problems, there are two possible solutions. 
First, the sample rate can be increased to produce more than two samples per pulse and 
more complicated detection routines can be employed. Another, simpler, solution is to 
set a truncation value. The first n sample pairs which exceed threshold are not used. The 
filter counter starts only after n detects. Since dropped pulses occur most frequently at 
the beginning and end of a pulse burst, the last n pairs in the pulse burst are also 
discarded. This method significantly reduces the number of dropped pulses. At present, 
the truncation value must be set by the operator. Ideally, this would be set by a CFAR 
(constant false-alarm rate) method [Ref. 4], The truncation value doubles as a flag that 
the burst is over. When n expected detects have been missed, processing of the data is 
suspended until just before the next, expected burst. 

It might seem that we are relying on a substantial amount of a priori knowledge 
about the signal to make this method work. Often there are specific known signals which 
are of interest, and accurate tailoring of the detection scheme to the emitter signal is an 
effective way to improve its performance. In addition, all that is really required is a good 
estimate of PRI. This comes from a single batch FFT process on the initial gain of a 
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emitter. In the future, the batch FFT could be used repetitively as a second check on the 
Kalman estimators. 

We are performing analog-to-digital conversion of the envelope detected signal at 
a rate of two to three times the reciprocal of the pulse width. However, other than simple 
threshold logic, we are not processing data at this rate. The algorithm actually samples 
the digital pulse train at only two to four times the PRF once it has been synchronized 
with the pulses. This method of PRI measurement appears promising with proper glitch 
rejection and dropped pulse logic. 

B. SCAN RATE MEASUREMENT 

Scan rate measurement stems from the PRI detection method. The First detect 
marks the beginning of a burst; the last marks the end. Processing these numbers allows 
calculation of the burst duration (and, therefore, the azimuthal bandwidth of the emitter). 
The distance from the start of one burst to the start of the next can be processed to give 
scan period. We are only looking for acquisition radars which have circular scans in 
azimuth. The scan rate is the number of 360° scans per unit time. Typical values for this 
would be 1 scan per second. 

As stated, the data must be sampled at a little more than twice once over the pulse 
width. This ensures at least one sample inside each pulse in the burst. A higher sampling 
frequency allows for more samples per pulse and would allow for pulse shape analysis. 
However, this is not the goal of this method, and the desire to minimize the amount of 
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data handling is strong. Minimizing the data required on an individual emitter, increases 
the number of possible emitters which the system could process. 

The time is also retained when the pulse starts. The next burst determines the initial 
estimate of scan rate. Scan rate is a constant which can be modelled as a digital system 
using the same system model as PRF: 



v *+i 



= x k ; 



= z k + 



( 13 ) 



where 

x k = scan rate estimates 

z k = scan rate observations, and 

v k = observation noise with variance r 0 . 

The development proceeds as before, and the final form of the filter is 

**♦1 = ** + y + t (z * " (14) 

Scan rate is assumed to be constant here. If the scan rate changes, it would be easy to 
impose limiting values for the change and prompt the operator when these are exceeded. 
The change could signal switch from scan to track mode as discussed earlier, a lost 
emitter contact, or the gain of a new emitter. 

Once the first estimate of scan rate has been made, the algorithm no longer needs 
to test every bit coming in. The program can use the initial estimate to narrow the range 
of data it processes in order to make less computation necessary to sort information from 
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the data. At this point the algorithm turns itself on at 90% of the expected scan time and 



begins processing again. 



C. TIME OF ARRIVAL 

During the PRF and scan rate measurement, it is possible to use an additional 
Kalman filter to produce a smoothed approximation to the burst envelope. The two sample 
detection scheme to identify and count pulses provides us with a means to begin to 
predict the shape of the pulse envelope. On initial detection of a pulse burst, the filters 
are initialized with discrete counters in order to allow processing on a digital computer. 
The observation vector consists of the two discrete samples which were used to detect the 
pulse. The magnitudes of the samples are summed and treated as a single estimate. The 
filter, since we are using a first order approximation, is developed in the same manner as 
before. This procedure of summing the two samples used to detect a pulse to form an 
observation z k is an ad hoc development. The final form is 



* k + 1 



= x k + 



V r o 

k+l/r. 






(15) 



where 

z k = individual samples in the pulse stream, 
r 0 = the noise variance, and 

x k = estimates of envelope shape. 

The counter k is only updated when observations occur, i.e., two consecutive samples 
exceed threshold. These estimates of envelope must be retained for the duration of a 
pulse burst to allow the TOA algorithm to calculate the energy centroid. The set of 
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estimates provides an approximation to the envelope which is integrated numerically in 
order to find the centroid of the shape. This centroid is what gives us the estimate of 
signal time of arrival. 

D. SNR AND THRESHOLD DETERMINATION 

To provide a measure of the signal-to-noise ratio is very important. If intelligence 
regarding the source strength of the emitting platform is known, then the signal strength 
can be used as a first approximation to the emitter range. Prior to the initial detect 
sequence, FFT methods are used to process the noise. The discrete Fourier transform is 
defined as [Ref. 6] 

N - 1 

X(k) = £ x{n)e- J(2nlN)nk . (16) 

n = 0 



The zero order coefficient is 

X(0) - E X(n). (17) 

n= 0 



Dividing by the number of samples gives us the mean: 



mean 



= 222 > = 



N 



N - 1 



E 4»> 

n-0 

N 



( 18 ) 



This method is convenient since the FFT of the data is also used to obtain initial 
estimates of the PRF. The mean along with the noise variance, can be used to set a 
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threshold value. A good threshold can be selected using detection theory concepts (i.e., 
CFAR detection): 

threshold = mean + a r 0 ; (19) 

where 

a = a floating value set by CFAR methods, and 
r 0 = the noise variance. 

The FFT is next used to calculate the power spectral density of the data: 

= x(*LxX(*r (20 ) 

N 



The total power in the data is 

P «m * E **> • <2I) 

k=0 

To estimate the SNR, the noise and signal powers are found separately. FFT's are taken 
in batches of 1024 for the noise only segments and 8192 for the signal segments. The 
measured and reported SNR is 



SNR = 10 



log 



signal 



8 xP, 



noise / 



( 22 ) 



These measurements and methods form the preprocessing in the area of extracting real 
time position information on emitting emitters and using Kalman optimal filtering 
techniques in order to provide emitter locations. 
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V. PERFORMANCE RUNS AND RESULTS 



This chapter shows graphical results of the methods described here. All the 
simulations were performed on a personal computer or a UNIX workstation. The real 
implementation of this technique would be such that, data enters as a stream of digitally 
sampled values from an analog to digital converter. The memory requirements are small, 
and only enough samples need be retained for the 8192 samples need to be retained in 
the batch FFT. To simulate this situation on a personal computer requires that a variable 
be loaded with the values for received data or that the computer manufacture the data "on 
the fly." Since all of the signal data is stored, the first situation can lead to memory 
problems during simulation which would not arise in a real implementation of this work. 
Even at moderate sample rates, the total amount of data used by the method quickly 
becomes large . The second method does not allow for displaying the data to reconstruct 
events which caused the algorithm to produce poor results. 

Due to the limitations described above, the emitter parameters were manipulated 
to improve the efficiency of the simulations. To permit more scans in a shorter period 
of time (fewer total data samples), the scan rates were set much higher than would be 
reasonable . The beamwidth of the emitter was similarly increased to provide the proper 
number of pulses per burst. Receiver velocity has been exaggerated to allow receiver 
movement to be seen during the short time span of the simulations. Therefore, values for 
parameters other than PRI will not be realistic. 
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A. EMITTER WITH ZERO NOISE 



To test the algorithm, a zero noise emitter was created. If the algorithm does not 
accurately predict the parameters of a noise-free pure signal, there is a fundamental error 
which must be corrected. Figure 9 shows the data used. There are five pulse bursts 
present. This data is the pulse train which will be seen by the first receiver. Data for the 
second receiver has the individual bursts shifted in time. The amount of shift is based 
on the emitter and receiver geometry. Factors such as Doppler shift due to the moving 
receiver are ignored. Data is manufactured so the precise time delays necessary for good 
angles of arrival exist in the data. If the algorithm finds the centroids properly, then the 
bearing lines to the emitter will be exact. 

The data in Figure 9 serves as the basis for all the tests presented here. For the 
case where noise is present, the data simply had zero mean, Gaussian noise added. The 
important parameters are: 

• Sample Frequency - 400 kHz 

• PRF - 9000 Hz 

• Scan Rate - 35 Scans/second 

• Beamwidth - 200° 

• Receiver velocity - 300,000 km/hr 

The test on the zero noise emitter was successful and plots of the results can be 
seen in Figures 10-12. Since there is no noise, the algorithm uses a threshold and a 
truncation value of one. There will not be any dropped pulses in this emitter. Figure 10 
shows the PRI history at both receivers. There is a bias in the Kalman filter of 
approximately 3 Hz. The source of the bias is, as yet, unknown. The final value of PRF 
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Figure 9. Emitter with Zero Noise 
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Figure 10. PRF Histories for the Zero Noise Emitter 
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Figure 11. AO A Estimates for Zero Noise Emitter 
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Figure 12. Closeup of AOA Estimates for Zero Noise Target 
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returned by the algorithm was 9003 Hz. This is better than the estimate from the batch 
FFT process. This estimate was 8984 Hz. This value served as the initial estimate for 
the Kalman filter. The reason for the FFT's error can be understood by looking at the 
resolution of the FFT: 

^ f _ Sample Frequency ^ , ^ 3 ) 

N 

where 

N = the number of samples in the batch FFT. 

The Af for the FFT of 8192 sample used is about 48 Hz. The Kalman filtering method 
provides superior resolution even with the current bias problem. The PRF estimate comes 
from looking for the maximum value in the FFT output. Some of the error from the FFT 
can be removed by averaging the results of FFT bins near the maximum. 

Scan rate histories are not shown. With only 4 to 6 pulse bursts in a typical 
simulation run, there is not much data to be seen. The probability of missing a burst is 
zero (for the simulation data) and the bursts are far enough apart in time, that the dropped 
pulse problem does not seriously affect the scan rate estimates. 

Figures 11 and 12 show the results of the angle of arrival (AOA) portion of the 
algorithm. Receiver positions over time are shown as circles in the lower portion of the 
graph. The lines show the estimated angles of arrival to the emitter for each set of 
receiver positions. For this simulation it is assumed that the receivers are a fixed distance 
apart and are traveling at a constant velocity. Figure 12 is a closeup view of the bearing 
lines in the vicinity of the emitter. Obviously, the estimated angles are sufficient to 
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predict the emitter bearing. Due to the narrow crossing angles of the bearing lines over 
time, position estimates of the emitter are not possible for these simulations. Even at the 
huge simulation velocity, the emitter bearings do not change much over time. Longer 
simulations using a different method generating emitter data will be required to produce 
sufficient numbers of emitter angles for tracking. 

B. EMITTER WITH ADDED NOISE 

Figure 13 shows the same emitter data as in Figure 9, but with added noise. Figure 
14 shows the PRI history using the program's a calculated threshold of 3.118 (calculated 
by Equation (19) with a=3) and a truncation value of 8. The effect of dropped pulses is 
quite clear in the results of this run. Pulses detected at the trailing edge of a burst reset 
the burst processing algorithm and an entire scan period becomes an estimate of PRI, this 
accounts for the large spikes in the graph for receiver 2. 

Figures 15 and 16 shows two more trials on the same file with adjusted threshold 
and truncation settings. Figure 15 again uses the threshold from Equation (19), but a 
truncation value of 12. Results for receiver 1 are improved, but some event disrupting 
the method has not been truncated for receiver 2. For the final run, the threshold was 
boosted manually to 3.5 (the calculated setting was 3.118) and truncation value was left 
at 12. Good performance can be seen at both receivers. At present, adjustment of these 
values is not automatically performed by the method. A study of adjusting these values 
based on a CFAR method would be useful. 
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Figure 17 gives the estimated bearing lines. Due to insufficient noise rejection by 
the algorithm, the results are poor. The estimates are only good enough to locate the 
quadrant of the emitter and not its actual position. More work is needed on the filtering 
and integration process which produce the estimated burst envelope and its centroid. 
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Figure 13. Emitter With Gaussian Noise 
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Figure 14. PRF History for Trial #1 
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Figure 15. PRF History for Trial #2 
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Figure 16. PRF History for Trial #3 
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Figure 17. AOA Estimates for Noisy Emitter 
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VI. SUMMARY 



A. CONCLUSIONS 

Further refinement of the method presented here and more work on classified data 
is in progress. The PRI estimates produced by the method are very accurate. These 
estimates are highly sensitive to dropped pulses. With the proper choice of threshold and 
truncation values, false or dropped detects were eliminated. Scan rate estimates were also 
good. The one significant drawback of the work done here is that it does not allow for 
easy application of itself to a radar or geometry that is different. The technique is 
precisely tailored to a particular type of emitter. However, having special algorithms can 
greatly improve the ability to detect specific emitters. 

B. RECOMMENDATIONS 

The author's suggestions for other research and points of interest occur sporadically 
throughout the text. Some final comments will be made here. The SNR of simulated 
data tested was high. This allowed testing and simulation to occur more easily. No 
research into the minimum SNR that can be analyzed by these techniques has been 
performed. Since simulations were performed on rather small data sets, the batch FFT 
processing occurs only on the first detected pulse. It is the author's contention that the 
FFT of every pulse burst is not necessary. More research should be performed to 
determine how often to reset the algorithm and update the noise measurements. 
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The algorithm presented currently has very poor resolution for the angle of arrival 
of the signal at low SNR. The method was developed with speed and minimum data 
storage as the primary goals. Improving the resolution could come by several means. 
Increasing the data rate will allow many more samples inside the burst and should 
increase the resolution at the cost of increased processing requirements. 
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APPENDIX A. ELINT RANGE ADVANTAGE 



Previously in this work, the fundamental range advantage of ELINT over radar was 
mentioned briefly, but never fully developed. This appendix will provide a more detailed 
approach to the range advantage of ELINT receivers. The equations and general flow of 
this derivation are from [Ref. 1], 

Consider the scenario of an emitting radar and an ELINT receiver attempting to 
intercept the radar signal. The signal power at a radar receiver is given by the basic radar 
equation [Ref. 1]: 



P t GtG r ^ 

(4k) 2 R^L t L r ’ 



(24) 



where 

S R = signal power at radar receiver (watts), 

P T = transmitted power (watts), 

G t = gain of radar transmit antenna, 

G r = gain of radar receive antenna, 

X = wavelength, 

a = emitter radar cross section (m), 

R r = radar losses from transmitter to antenna, 

L t = radar losses from the transmitter to the antenna, and 
L r = radar losses from antenna to receiver. 

The ELINT receiver, however, does not rely on the reflected energy from an emitter and 
so we use the one-way or "beacon" radar equation [Ref. 1]: 



„ P t G te G e^ . 
(47t) 2 R^L T L E ’ 



( 25 ) 
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where 



S E = signal level at the ELINT receiver (watts), 

^te = gain of radar transmit antenna in the direction of the ELINT receiver, 

G e = gain of the ELINT antenna, 

R e = range from radar to ELINT receiver, 

L e = loss from ELINT antenna to receiver, and 

all other symbols are as they appear in Equation (24) above. 

Suppose that the signal level required at the ELINT receiver is a factor of 8 times the 
signal level required at the radar to detect the emitter: 

S E = 65*. (26) 

Suppose, also, that the range, R E , at which the ELINT received power level is 5S R is a 
multiple of a times that range at which the radar received power level is S R : 

R e = a R r . (27) 

Substituting Equations (26) and (27) into Equation (25) above, gives [Ref. 1]: 



^ _ Pt-Gj^G^A, 

b(4n) 2 R 2 R a 2 L T L E 

By substituting Equation (28) into Equation (24) and solving for a we obtain 



(28) 
[Ref. 1] 



i = — = R 4 7t G TE G E L r 

R r R b a G t G r L e 

Suppose, for convenience, that the losses in the radar and ELINT receivers are equal. 
Suppose also that the radar’s transmit and receive antenna gains are equal and that the 
ELINT receiver is in the side lobe of the radar such that G te 1. Finally, suppose that 
the ELINT receiver antenna is omnidirectional (G E =1). Then Equation (29) becomes 



(29) 
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a 



2 



(30) 



Rr 4 k 
G r i 6 



The factor (l/5) 1/: is the advantage enjoyed by the radar due to its more sensitive receiver. 
(The radar receiver can be matched to its signal and can use multiple pulse integration.) 
The ELINT receiver generally has a wider noised bandwidth than the radar receiver and 
operates on a single pulse basis. Therefore the minimum that 5 can be is one and it is 
typically as large as 100 or 1000. The radar antenna gains are typically in the 30 to 40 
dB range. Figure 18 shows Equation (30) plotted for several values of 5 and G R . Figure 
18 also shows the main beam intercept case, in which G TE = G R for this case, the ELINT 
range advantage is given by 



a = R b 



4k 

6 O Gr 



• 2 
2 



(31) 



Let us assume that our ELINT receiver is in the side lobes of a radar having a 40 dB gain 
antenna, which can detect a 1 nr emitter at 100 km. If the ELINT receiver is 20dB less 
sensitive than the radar receiver (8 = 100) then the ELINT range advantage from Equation 
(30) is 



a 



10 5 m 


471 


10* 


100x(l m 2 ) 



2 _ 



« 3.5. 



(32) 



This is a substantial advantage over the radar. If we move the receiver into the main 
beam then, for the same radar and ELINT receiver, Equation (31) gives the new ELINT 
range advantage: 
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= 350. 



(33) 



a = XQr’m 



4n 



■_i 

2 



100 x(l m 2 )x 10 4 



This means that we should be able to detect the radar at a range of 35,000 kilometers. 
This is, of course, only a theoretical detection range and does not account for such effects 
as line-of-sight restrictions. 




Figure 18. ELINT Range Advantage 
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APPENDIX B. MATLAB SOURCE CODE 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



% % 

% WORLD.M % 

% Written by Gary Sipe % 

% % 

% % 

% This program is a menu-driven utility used to write the geometry files % 

% used by MAKETGT.M to construct target signal files. It is menu driven % 

% and its operation is self explantory. The filename given when requested % 

% by the program must no more than 5 characters in length. The data % 

% output by WORLD will be stored in the file FNAME-G.MAT where % 

% "FNAME" is the filename you specified at the prompt. % 

% % 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



clear; 

clg; 



% Initialize the default values of the system variables. 



samp_freq=400000; 

prf=9000; 

pw=5e-6; 

beam_w=7.5; 

scan_rate= 1 ; 

duration=0.05; 

xl=-25; 
y 1=0; 

x2=25; 

y2=0; 

xt=50; 

yt=100; 

v=[300 330]; 

menut='Setup Options'; 



%in Hz 
%in Hz 
%in seconds 
%in degrees 
%in sweeps/second 
%in seconds 

%in km 
%in km 

%in km 
%in km 

%in km 
%in km 

%in [km/hr bearing] 
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menu 1 -Change Radar Parameters'; 
menu2='Change Receiver Parameters'; 
menu3- Change Target Parameters'; 
menu4-Load Old File'; 
menu5='All Done'; 

done=0; 
format compact 

% Enter the main program loop. The loop does nothing more than display current 
% settings and prompt the user to make changes. Resulting changes may be 
% saved to a file or discarded. 

while done==0, 
clc; 

disp(['PRF = ' int2str(prf) ' Hz']); 

disp(['Pulse Width = ' int2str(pw*le6) ' microseconds']); 
disp(['Beam Width = 1 num2str(beam_w) ' degrees']); 
disp(['Scan Rate = ' num2str(scan_rate) ' sweeps/second']); 
disp(['Duration of Problem = 1 num2str(duration) ’ seconds']); 
disp("); 

disp('Receiver and Target Coordinates are in kilometers'); 

disp("); 

disp(['Sample Frequency = ' int2str(samp_freq) ' Hz']); 
disp(['Receiver 1 initial position = (' num2str(xl) num2str(yl) ')']); 

disp(['Receiver 2 initial position = (' num2str(x2) num2str(y2) ')']); 

disp(['Receiver velocity = ' num2str(v(l)) ' km/hr']); 
disp(['Receiver direction = ' num2str(v(2)) ' degrees true']); 
disp(['Target initial position = (' num2str(xt) num2str(yt) ')']); 

choice=menu(menut,menul,menu2,menu3,menu4,menu5); 
if choice==l, 
clc; 

new=input(['Input the New PRF (<ENTER>- int2str(prf) ') ']); 
if -isempty(new), 
prf=new; 
end; 

clear new 

new=input(['Input the New Pulse Width (<ENTER>- int2str(pw*le6) ') ']); 
if ~isempty(new), 
pw=new/le6; 
end; 

clear new 

new=input(['Input the New Beam Width (<ENTER>- num2str(beam_w) ') ']); 
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if ~isempty(new), 
beam_w=new; 
end; 

clear new 

new=input(['Input the New Scan Rate (<ENTER>=' num2str(scan rate) ') ']); 
if -isempty(new), 
scan_rate=new; 
end; 

clear new 

new=input(['Input the New Duration of Problem (<ENTER>- num2str(duration) ') 

); 

if -isempty(new), 
duration=new; 
end; 

clear new 
elseif choice==2, 
clc; 

new=input(['Input the New Sample Freq (<ENTER>- int2str(samp_freq) ') ']); 
if ~isempty(new), 
samp_freq=new; 
end; 

clear new 

new=input(['Input the xl (<ENTER>- num2str(xl) ') ']); 
if ~isempty(new), 
xl=new; 
end; 

clear new 

new=input(['Input the yl (<ENTER>- num2str(yl) ') ']); 
if -isempty(new), 
y l=new; 
end; 

clear new 

new=input(['Input the x2 (<ENTER>- num2str(x2) ') ']); 
if ~isempty(new), 
x2=new; 
end; 

clear new 

new=input(['Input the yl (<ENTER>- num2str(y2) ') ']); 
if -isempty(new), 
y2=new; 
end; 

clear new 

new=input(['Input the New Velocity (<ENTER>- num2str(v(l)) ') ']); 
if -isempty(new). 
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v(l)=new; 

end; 

clear new 

new=input(['Input the New Direction (<ENTER>- num2str(v(2)) ') ']); 
if -isempty(new), 
v(2)=new; 
end; 

clear new 
elseif choice==3, 
clc; 

new=input(['Input the xt (<ENTER>- num2str(xt) ') ']); 
if -isempty(new), 
xt=new; 
end; 

clear new 

new=input(['Input the yt (<ENTER>- num2str(yt) ') ']); 
if -isempty(new), 
yt=new; 
end; 

clear new 
elseif choice==4, 

fname_in=input('Input the problem name (5 chars/no extension) ','s'); 
fname_in=[fname_in *-g*]; 
eval(['load 1 fname_in]); 
elseif choice==5, 
done=l; 
end; 
end; 

done=0; 

ans=input('Save Changes to a Problem File? (y/n) ','s'); 
if ans=='y', 
while done==0, 

fname_out=input('Input the problem name (5 chars/no extension) ','s'); 
fname_out=[fname_out '-g']; 
old=exist([fname_out '.mat']); 
if old==2, 

ans=input('File exists — delete it? (y/n) ','s'); 
if ans=='y', 

clear menul menu2 menu3 menu4 menu5 done ans choice menut 
eval(['save ' fname_out]); 
done=l; 
end; 
else, 
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clear menul menu2 menu3 menu4 menu5 done ans choice menut 
eval(['save ' fname out]); 
done=l; 
end; 
end; 
end; 

format loose 



54 




% 

% MAKETGT.M 
% Written by Gary Sipe 



% 

% 

% 

% 



% This program reads the geometry and radar data provided by the '-G' files % 
% created with WORLD. M. Filenames are again only 5 letters in length. % 
% The signal file created by MAKETGT will be saved in FNAME-S.MAT. % 



% Reset the environment. 

clear 

clg 

rand('normal'); 

% Get the problem file of interest, 
clc; 

fname=input('Input the problem name (5 chars/no extension) 
eval(['load ' fname '-g']); 

ls=2.997925e5; %(kilometers/sec) 

total_l=duration*samp_freq; 

delay=5000; 

toa_al=delay; 

% Prompt regarding noise addition and initialize the data vectors accordingly. 
ans=input('Do you want to add noise? ’/s'); 
if ans=='y', 
data 1 =rand( 1 ,total_l) ; 
else, 

data 1 =zeros( 1 ,total_l) ; 
end; 

data2=datal; 

d=length(datal); 

disp('Working'); 

% Calculate and build a burst envelope which will then be displaced in time in the 

% vectors DATA1 and DATA2 to build the final data vectors. 

nextscan=floor(samp_freq/scan_rate); 

nexttime=nextscan/samp_freq; 

numscans=ceil((total_l-delay)/nextscan); 

target_l=floor((beam_w/360/scan_rate)*samp_freq); 



% 



% 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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envelope=sin(linspace(0,pi,target_l)); 
time=linspace(0,beam_w/360/scan_rate,target_l); 
mod=(square(2*pi*prf*time,pw*100*prf)+ones(time)).*0.5; 
target=envelope.*mod.* 1 0; 

% Build data vectors using the burst envelope created above, 
for c=l:numscans, 

xd=sqrt((x2(c)-xt) A 2+(y2(c)-yt) A 2)-sqrt((xl(c)-xt) A 2+(yl(c)-yt) A 2); 
xdn=ro un d(xd * sam p_freq/l s) ; 
if xdn+target_l+toa_al(c)<=d, 

datal(l,toa_al(c):toa_al(c)+target_l-l)=datal(l,toa_al(c):toa_al(c)+target_l-l)+target; 

toa_a2(c)=toa_al(c)+xdn; 

data2(l,toa_a2(c):toa_a2(c)+target_l-l)=data2(l,toa_a2(c):toa_a2(c)+target_l-l)+target; 
toa_a 1 (c+ 1 )=toa_a 1 (c)+nextscan; 
end; 

if c~=numscans, 

xl(c+l)=xl(c)+v(l)*nexttime*sin(v(2)*pi/180)/3600; 
x2(c+l)=x2(c)+v(l)*nexttime*sin(v(2)*pi/180)/3600; 
y l(c+l)=y l(c)+v(l)*nexttime*cos(v(2)*pi/l 80)/3600; 
y2(c+l)=yl(c+l); 
end; 
end; 

clc; 

% Prompt user for filename and save the necessary data. 

% Extraneous variables are cleared prior to saving. 
ans=input('Do you want to save this target? (y/n) ','s'); 
if ans == y, 

clear c ans mod target target_l target_t toa_a toa_a2(c) nexttime 
clear numscans con Is d delay xdn totalj duration 
clear time nextscan envelope 
eval(['save ' fname ’-s']); 
end; 
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% 

% MAIN.M 
% Written by Gary Sipe 

% 



% 

% 

% 

% 



% This program is the central processing loop of the algorithm. It calls the % 
% subroutine PROCESS to perform the processing of individual pulse bursts. % 
% it will retrieve a signal file ('fname-s.mat') and output all parameters to an % 



% Retrieve the problem file of interest, 
clc; 

fname=input('Input the problem name, 
loadstr=['load ' fname '-s']; 
eval(loadstr); 

time=linspace(0,length(datal)/samp_freq,length(datal)); 

new=[]; 

% The next series of program blocks prompts the user for various options. 
ans=input('Do you want to view the target? (y/n) ','s'); 
if ans == 'y', 

plot(time,datal),grid,xlabel('Time -- seconds'),ylabel('Magnitude'), 
title(['Plot of Data for Target File -- ’ fname]) 



ans=input('Do you want to view the intermediate results? (y/n) ','s'); 
if ans == *y'» 
show=l; 
else, 
show=0; 
end; 

ans=input('Do you want keyboard access between results? (y/n) ','s'); 
if ans == 'y', 
key=l; 
else, 
key=0; 



% analysis file ('fname-a.mat'). 



% 




% Reset the environment, 
clear 



clg 

format compact 



pause; 

end; 



57 



end; 



clc; 

disp('Beginning Analysis — 
disp("); 

[N _mean,N_var,N_first,N_totp wr]=my fft(data 1(1:1 024),samp_freq); 

[S _mean,S_var,S_first,S_totpwr]=myfft(data 1(500 1:131 92),samp_freq); 

a=3; 

pri=S_first; 
if N_totpwr 0, 

SNR=10*logl0(S_totpwr/8/N_totpwr); 
th resh =N_m ean+a* N_var; 
else, 

thresh=l; 

SNR='undefined'; 

end; 

disp(['Auto threshold is — ' num2str(thresh) ' <ENTER> to accept']); 
new=input('Or type in your own threshold now'); 
if -isempty(new), 
thresh=new; 
end; 

for sens=l:2; 
if sens==l, 

disp('Receiver 1 Simulation'); 
data=datal; 
else, 

disp('Receiver 2 Simulation'); 
data=data2; 
pri=S_first; 
end; 

disp('In Loop 1'); 
for cl = l:l:length(data), 
if data(cl)>thresh, 
if data(cl+l)>thresh, 
disp('detect'); 
scan=l ; 

[pri,xhat,start,done]=process(data,cl,samp_freq,pri,thresh, show, key); 
[toa(scan)]=centroid(xhat,samp_freq); 
break; 
end; 
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end; 

end; 

disp('in Loop 2') 
done=done+l 000; 
lastscan=start; 
for c2=done:l:length(data), 
if data(c2)>thresh, 
if data(c2+l)>thresh, 
disp('detect'); 
start=c2; 
scan=2; 

[pri,xhat, start, done]=process(data,c2,samp_freq,pri,thresh,show,key); 
[toa(scan)]=centroid(xhat,samp_freq); 
scan_rate( 1 )=samp_freq/(start-lastscan); 
lastscan=start; 
break; 
end; 
end; 
end; 



c3=done+1000; 
finished=0; 
disp('In Loop 3') 

while finished==0, 

if c3>=length(data); 

break; 

else; 

c3=c3+l; 

end; 

if data(c3)>thresh, 
if data(c3+l)>thresh, 
disp('detect'); 

[pri,xhat,start,done]=process(data,c3,samp_freq,pri,thresh,show,key); 

scan_obs=samp_freq/(start-lastscan); 

scan_rate(scan)=scan_rate(scan-l)+(scan_obs-scan_rate(scan-l))/scan; 

scan=scan+l; 

lastscan=start; 

[toa(scan)]=centroid(xhat,samp_freq); 

c3=done+1000; 
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end; 

end; 

end; 



if sens==l, 
toal=toa; 
pril=pri; 

scan_rate 1 =scan_rate ; 
scan_rate=0; 
else, 

toa2=toa; 

pri2=pri; 

scan_rate2=scan_rate; 

end; 

end; 

clear data datal data2 time 
eval(['save ' fname '-a']); 
format loose 
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%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% % 

% MYFFT.M % 

% Written by Gary Sipe % 

% % 

% This function is called by MAIN to obtain spectrum information during % 

% the first detection of an emitter. It produces the average value, the % 

% variance, the first spectral line, and the total power. % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



function [avg,var,first,totpwr]=myfft(data,samp_freq); 
%function [avg,var,first,totpwr]=myfft(data,samp_freq); 

l=length(data); 

y=fft(data); 

avg=y(l)/l; 

spec_den=(y.*conj(y)/l); 

totpwr=sum(spec_den); 

var=cov(data); 

delf=400000/l; 

[dummy, first]=max(y(2: 250)); 
first=(first-l)*delf; 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% % 

% PROCESS. M % 



% Written by Gary Sipe % 

% % 



% This function is called by MAIN to handle burst processing. It returns the % 
% estimates of pri, the burst envelope, the burst start time, and the burst % 
% ending time. This function also handles all truncation logic. % 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



function [pri,xhat, start, done]=process(data,start,samp_f,pri,thresh, show, key); 

%function [pri,xhat, start, done]=process(data, start, samp_f,pri,thresh,show, key); 

% The variables 'show' and 'key' are flags used to allow displays during the processing 
% cycle, 
if nargin==5, 
show=0; 
else 

show=show; 

end; 



if nargin==6, 
key=0; 
else, 

key=key; 

end; 

% Set all flags and inital values. 

tau=l/samp_f; 

broken=0; 

xhatold=0; 

trunc=12; 

cold=zeros(l,trunc); 

a=0; 

p=length(pri); 

k=0; 

pri_old=pri(p); 

nuls=0; 

disp('In process'); 

for c=start:l:length(data), 
nuls=nuls+l; 
if data(c)>thresh, 
if data(c+l)>thresh, 
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nuls=0; 

if cold(trunc)~=0, 

% Truncate the first few pulses based on the 'trunc' setting, 
if a<=trunc, 
a=a+ 1 ; 
end; 

if a==trunc, 
a=a+l; 
start=c; 
end; 

% Beginning of burst processing routine, 
if a> trunc, 

% PRI logic and filter. 
p=p+l; 

pri_hat=(samp_f/(c-cold(trunc))); 

pri(p)=pri_old+(pri_hat-pri_old)/(p+l); 

pri_old=pri(p); 

% Envelope logic and filter. 
k=k+l; 

xhat(k)=xhatold+(l/k)*((data(c)+data(c+l))-xhatold); 

xhatold=xhat(k); 

end; 

end; 

% Simulate a set of shift registers to be used when truncating the end of the burst. 
cold=shiftl(cold,c); 
end; 
end; 

% This is the end of burst logic, 
if a>trunc, 

if nuls>trunc*samp_f/pri(p); 
done=cold(l); 

pri=pri( 1 :(length(pri)-trunc)); 
xhat=xhat( 1 :(length(xhat)-trunc)); 
xhat(length(xhat)+l )=start; 
xhat(length(xhat)+l)=cold(l); 
broken=l; 
break; 
end; 
end; 
end; 
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% Exit logic, 
if broken==0, 
done=length(data); 
pri=pri( 1 .(length(pri)-trunc)); 
xhat=xhat(l:(length(xhat)-trunc)); 
xhat(length(xhat)+l )=start; 
xhat(length(xhat)+ 1 )=cold( 1 ); 
end; 

% Display logic, 
if show~=0, 

subplot(21 l),plot(pri),grid 
subplot(212),plot(xhat(l:length(xhat)-2)),grid 
pause; 
clg; 
end; 

% Keyboard utility logic, 
if key==l, 
keyboard; 
end; 
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% This function is called by MAIN to numerically integrate the estimated 
% burst envelope. 



% 

% CENTROID. M 
% Written by Gary Sipe 
% 



% 

% 

% 

% 

% 

% 




function [toa]=centroid(xhat,samp_freq); 

%function [toa]=centroid(xhat); 

l=length(xhat); 

minx=xhat(l-l); 

maxx=xhat(l); 

xhat=xhat(l:l-2); 

M=0; 

Mx=0; 

for count=l:l-3, 

M=M+0.5*(xhat(count)+xhat(count+l)); 

Mx=Mx+0.5*(count+0.5)*(xhat(count)+xhat(count+l)); 

end; 

x_bar=Mx/M; 

toa=minx+x_bar*(maxx-minx)/(l-2); 

toa=toa/samp_freq; 
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%%%%%%% 0 /o 0 /o 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%% 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%%% 
% % 

% SHIFT.M 

% Written by Gary Sipe 

0 / 0 / 

/o /o 

% Simulates a block of shift registers. % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



% 

% 



function [data]=shift(data,new); 

l=length(data); 

data(l:l-l)=data(2:l); 

data(l)=new; 

return; 
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%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%% 0 /o 0 /o%% 
% % 

% BIRDSEYE. M % 

% Written by Gary Sipe % 

% % 

% This program file is used to provide the graphical displays of estimated % 

% lines of position over time. It reads the analysis files (fname-a.mat). % 

0 /o 0 /o% 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%% 0 /o%% 0 /o 0 /o 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%%% 



clear 

clg 

ls=2.997925e5; % Speed of light (km/sec) 

fnameHnputOInput the problem file name 
eval(['load ' fname '-a']); 

axis('square') 

axis([-50 150 0 350]); 

plot(xt,yt,'o') 

hold on; 

xoff=0; 

yoff=0; 

for c=l:length(toal), 
dt=toal (c)-toa2(c); 

dtmax=sqrt((xl(c)-x2(c)) A 2+(y l(c)-y2(c)) A 2)/ls; 

cnt=sqrt((xl(c)-x2(c)) A 2+(yl(c)-y2(c)) A 2)/2; 

a=cnt*dt/dtmax; 

b=sqrt(cnt A 2-a A 2); 

xoff=(xl(c)+x2(c))/2; 

yoff=(y l(c)+y2(c))/2; 

x=linspace(0,500,1000); 

LOP=(x)*b/a+yoff; 
plot(xl(c),y l(c),'rx'),grid 
plot(x2(c) > y2(c),'rx , ),grid 
plot(x+xoff,LOP,'g-') 
pause(3); 
end; 

xlabel('X in km.’), 
ylabel('Y in km.'), 
hold off; 
axis('normar); 



67 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



% % 

% EXAMINE. M % 

% Written by Gary Sipe % 

% % 

% This utility is used to interactively display and save the contents of very % 

% variables. Screen buttons are drawn on the MATLAB display and the % 

% contents can be scrolled. Many other options are also available such as % 

% plot labels, changing axes, changing starting point and block sizes of the % 

% display, and saving the current screen to a meta file. This is a utility % 

% function used during debugging and no documentation is provided. % 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



function examine(time,data); 

clg 

clc 



if nargin==l, 
data=time; 

time=linspace(l,length(data),length(data)); 

end; 

maxval=max(data); 

minval=min(data); 

l=length(data); 

disply=l; 

hardcopy=0; 

getname=T; 

reverse=0; 

forward=l, 

tstring-'; 

xstring- 

ystring="; 

menutitle-Labels Menus'; 
menuI-Change Title'; 
menu2='Change Xlabel'; 
menu3- Change Ylabel'; 
menu4='Place a Text Marker'; 
menu5='Preview Your Plot'; 
menu6- Save Plot to a Meta File'; 
menu7='Back to Plot'; 
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clc 

ans=input('Will you want to save plots to a meta file? (y/n) ','s'); 
if ans=='n', 
nomet=l ; 
else, 

while getname==l; 
nomet=0; 
clc 

disp(’Input the filename for meta files (8 chars no extension) '); 
metaname=input('8 chars no extension (n=none) ','s'); 
metaname=[metaname '.met']; 
ans=exist(metaname); 
if ans==2; 

ans=input('File exists -- delete it? (y/n) ','s'); 
if ans=-y' 

eval(['!del ' metaname]); 
getname=0; 

eval(['meta ' metaname]); 
end; 
else, 

getname=0; 

eval(['meta ' metaname]); 
end; 
end; 
end; 

a=l; 

block=100; 

butdx=[0.15 0.215 0.215 0.15 0.15]; 
butlx=[0.26 0.35 0.35 0.26 0.26]; 
butax=[0.35 0.44 0.44 0.35 0.35]; 
buthcx=[0.44 0.53 0.53 0.44 0.44]; 
butsx=[0.575 0.665 0.665 0.575 0.575]; 
butbx=[0.665 0.755 0.755 0.665 0.665]; 
butrx=[0.80 0.85 0.85 0.80 0.80]; 
butfx=[0.85 0.90 0.90 0.85 0.85]; 
buty=[0.938 0.938 0.888 0.888 .938]; 
axis([time(a) time(a+block) minval maxval]) 



69 



% The remainder of the program is the main loop. Tracing the operation of this loop 
% and the logic required for the screen buttons is not recommended. 

while disply==l; 
if hardcopy— 0, 

axis([time(a) time(a+block) minval maxval]) 
displot(time(a:a+block),data(a:a+block)),grid 
title(tstring),xlabel(xstring),ylabel(ystring); 
polyline(butdx,buty,'sc') 
text(0. 1 55,0. 898, 'Done', 'sc') 
polyline(butlx,buty,'sc') 
text(0.265,0.898,'Labels',’sc') 
polyline(butax,buty,'sc') 
text(0. 37,0. 898, 'Axis', 'sc') 
if nomet==0, 
polyline(buthcx,buty,'sc') 
text(0. 46,0. 898, 'Meta'.'sc') 
end; 

polyline(butsx,buty,'sc') 
text(0.5 85,0.898, 'Start','sc') 
polyline(butbx,buty,'sc') 
text(0.675,0.898,'Block','sc') 
if reverse==l, 
polyline(butrx,buty,'sc') 
text(0. 805,0. 898, 'Rev','sc') 
end; 

if forward— 1, 
polyline(butfx,buty,'sc') 
text(0.855,0.898,'Fwd','sc') 
end; 

[x,y]=ginput(l,'sc'); 

axis; 

if x>0.80&x<0.85&y>0.888&y<0.938, 
if reverse==l, 
if a-block<=l; 
a=l; 

reverse=0; 

forward=l; 

else; 

a=a-block; 

reverse=l; 

forward=l; 
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end; 

end; 

elseif x>0.85&x<0.90&y>0.888&y<0.938, 
if forward==l, 
if a+block>=l-block, 
a=l-block; 
reverse=l; 
forward=0; 
else, 

a=a+block; 

reverse=l; 

forward=l; 

end; 

end; 

elseif x>0. 1 5&x<0.21 5&y>0.888&y<0.938, 
disply=0; 

elseif x>0.26&x<0.35&y>0.888&y<0.938, 
loop=l; 

while loop==l; 
clc 

disp(['Current title is: ' tstring]) 
disp(['Current Xlabel is: ' xstring]) 
disp(['Current Ylabel is: ' ystring]) 
disp(") 
disp(") 

choice=menu(menutitle,menul,menu2,menu3,menu4,menu5,menu6,menu7); 
if choice==l, 

tstring=input('Input the new title 
elseif choice==2, 

xstring=input('Input the new Xlabel ','s'); 
elseif choice==3, 

ystring=input('Input the new Ylabel ','s'); 
elseif choice==4, 

tmstring=input('Input the string for the text marker ','s'); 

axis([time(a) time(a+block) 0 maxval]); 

displot(time(a:a+block),data(a:a+block)),grid 

title(tstring),ylabel(ystring),xlabel(xstring) 

gtext(tmstring), 

pause(2); 

axis; 

elseif choice==5, 

axis([time(a) time(a+block) 0 maxval]); 

displot(time(a:a+block),data(a:a+block)),grid 

title(tstring),ylabel(ystring),xlabel(xstring) 
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pause; 

axis; 

elseif choice==6, 
meta; 

elseif choice==7, 
loop=0; 
end; 
end; 

elseif x>0.35&x<0.44&y>0.888&y<0.938, 
clc 

disp('Values for the x-axis are scaled by the plot window ') 
disp('based on the starting point and block size. You may only ') 
disp('scale the y-axis of the plot.') 
disp(") 
disp(") 
new =[]; 

new=input(['Input the new y-min (<ENTER>=' num2str(minval) ') ']); 
if -isempty(new), 
minval=new; 
end; 
new=[]; 

new=input(['Input the new y-max (<ENTER>- num2str(maxval) ') ']); 
if -isempty(new), 
maxval=new; 
end; 

elseif x>0.575&x<0.665&y>0.888&y<0.938, 
new=[]; 
clc 

new=input(['Input the new starting point (<ENTER>- int2str(a) ') ']); 
if -isempty(new), 
a=new; 
reverse=l; 
forward=l; 
end; 

if a<=l, 
a-l; 

reverse=0; 

forward=l; 

elseif a+block>=l-block, 
a=l-block; 
reverse=l; 
forward=0; 
end; 

elseif x>0.665&x<0.755&y>0.888&y<0.938, 
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nevv=[]; 

clc 

new=input(['Input the new block size (<ENTER>- int2str(block) ') ']); 
if ~isempty(new), 
if new>=l, 
block=l-l; 
a=l; 

forward=0; 
reverse=0; 
elseif new<l, 
block=l; 
else, 

block=new; 

end 

end; 

if a+block>=l-block, 
a=l-block; 
reverse=l; 
forward=0; 
end; 

elseif nomet==0, 

if x>0.44&x<0.53&y<9.38&y>0.888, 
if hardcopy — 1; 

hardcopy = 0; 
else; 

hardcopy=l; 

end; 

end; 

end; 

elseif hardcopy— 1, 

displot(time(a:a+block),data(a:a+block)),grid 

title(tstring),ylabel(ystring),xlabel(xstring) 

meta; 

text(0.40,0.7,'Press Any Key to Continue’,'sc') 
pause; 

hardcopy=0, 

axis; 

end; 

end; 

disp("); 

clc; 
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